function [waveformShape,RFE,RFEbyFreq,holdoutRFE,holdoutRFEbyFreq]=multiSpikeTrainToLFP(lfp,spikes,parameters)
spikes=bsxfun(@minus,spikes,mean(spikes));
%% Extract parameters
lowFreq=parameters.lowFilterFrequency;
highFreq=parameters.highFilterFrequency;
KernelWidth=parameters.dicionaryTimeSpan;
frequencies=parameters.testingFrequencies;
frequencyBandwidth=parameters.testingFrequencyBandwidth;
samprate=parameters.samplingRate;
numCells=size(spikes,2);
pen=parameters.ridgeRegression;
T=numel(lfp);
%% Filter the initial LFP data
if highFreq>50
    filter60HzSignal=true;
else
    filter60HzSignal=false;
end
lfpFilt=filterLFP(lfp',samprate,lowFreq,highFreq,filter60HzSignal)';
%% Calculate Convolutional Kernels
RFE=zeros(numCells,1); % Reduction in Fractional Error
convKernels=zeros(2*KernelWidth+1,1); % Convolutional Kernels
[RFE,convKernels]=timeInvariantMultiCell(lfpFilt,spikes,KernelWidth,pen);
waveformShape=convKernels;
%% Calculate LFP predictions from Single Spike Trains
lfpPred=zeros(numel(lfpFilt),1);
for n=1:numCells
    tmp=conv(spikes(:,n),convKernels(:,n),'full');
    lfpPred=lfpPred+tmp(KernelWidth+1:end-KernelWidth);
end
%% Calculate RFE by Frequency Bins
[RFEbyFreq]=calcRFEbyFrequency(lfpFilt,lfpPred,...
    samprate,frequencies,frequencyBandwidth);

%% cross-validation
cv=parameters.numberOfCrossValidationRuns;
if cv==0
    holdoutRFE=[];
else
    Tbin=floor(T/cv);
    lfpPredcv=zeros(Tbin*cv,1);
    for c=1:cv
    fprintf('CV #%d\n',c);
    ndx=[1:(c-1)*Tbin,Tbin*c+1:Tbin*cv];
    lfpcv=lfpFilt(ndx);
    spikescv=spikes(ndx,:);
    [~,convKernelscv]=timeInvariantMultiCell(lfpcv,spikescv,KernelWidth,pen);
    %%
    mi=(c-1)*Tbin+1-KernelWidth;
    ma=Tbin*c+KernelWidth;
    predndx2=mi:ma;
    predndx=[(c-1)*Tbin+1:Tbin*c];
    lfpPredtmp=zeros(Tbin,1);
    tmp2=zeros(Tbin+2*KernelWidth,1);
    for n=1:numCells
        tmp=conv(spikes(predndx,n),convKernelscv(:,n),'full');
        tmp2=tmp+tmp2;
%         lfpPredtmp=lfpPredtmp+tmp(KernelWidth+1:end-KernelWidth);
    end
    use=(predndx2>0) & (predndx2<Tbin*cv);
    lfpPredcv(predndx2(use))=lfpPredcv(predndx2(use))+tmp2(use);
    
    end
    holdoutRFE=1-sum((lfpFilt(1:Tbin*cv)-lfpPredcv).^2)./sum((lfpFilt(1:Tbin*cv)).^2);
    [holdoutRFEbyFreq]=calcRFEbyFrequency(lfpFilt(1:Tbin*cv),lfpPredcv,...
    samprate,frequencies,frequencyBandwidth);
    1;
end
return
%%
ndx=5e3:5e3+100;
plot(ndx,lfpFilt(ndx),'k-',ndx,lfpPred(ndx),'b-.',ndx,lfpPredcv(ndx),'r--','LineWidth',2)
%%
ndx=5e3:5e3+100;
plot(ndx,lfpFilt(ndx),'k--',ndx,lfpFilt(ndx)-lfpPred(ndx),'r-','LineWidth',2)
%%
[acf,lags,bounds] = crosscorr(lfpPred(1:Tbin*cv),lfpPredcv,500)
plot(lags,acf)











